Isomorphs in model molecular liquids 
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Isomorphs are curves in the phase diagram along which a number of static and dynamic quan- 
tities are invariant in reduced units. A liquid has good isomorphs if and only if it is strongly 
correlating, i.e., the equilibrium virial/potential energy fluctuations are more than 90% correlated 
in the NVT ensemble. This paper generalizes isomorphs to liquids composed of rigid molecules 
and study the isomorphs of two systems of small rigid molecules, the asymmetric dumbbell model 
and the Lewis— Wahnstrom OTP model. In particular, for both systems we find that the isochoric 
heat capacity, the excess entropy, the reduced molecular center-of-mass self part of the intermedi- 
ate scattering function, the reduced molecular center-of-mass radial distribution function to a good 
approximation are invariant along an isomorph. In agreement with theory, we also find that an 
instantaneous change of temperature and density from an equilibrated state point to another iso- 
morphic state point leads to no relaxation. The isomorphs of the Lewis— Wahnstrom OTP model 
were found to be more approximative than those of the asymmetric dumbbell model, which is 
consistent with the OTP model being less strongly correlating. For both models we find "master 
isomorphs", i.e., isomorphs have identical shape in the virial/potential energy phase diagram. 
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I. INTRODUCTION 

For supercooled liquids near the glass transition chang- 
ing slightly the density p or temperature T the struc- 
tural relaxation time may change several orders of 
magnitude. In the study of these liquids^"— it is found 
that Ta does not change when p'^ /T is kept constant, 
where 7 is a fixed material-specific exponent. This phe- 
nomenon is called density scaling and has been estab- 
lished for many liquids, excluding associative liquids such 
as water~. A related observation is isochronal superpo- 
sition^^'^, i.e., that supercooled state points with iden- 
tical Tq, have the same dielectric spectrum. A different 
and at first sight unrelated concept is Rosenf eld's excess 
entropy scaling^. In this procedure a relation is estab- 
lished between hard-to-predict dynamic properties and 
easier-to-predict thermodynamic quantities, here the ex- 
cess entropy via a scaling of the dynamics to so-called re- 
duced units. Initially this was observed for model atomic 
liquids^ii, but later it was extended to model molecular 
liquids^"— and experimental liquidsii. 

In a recent series of paper&i^r— a new class of liquids 
was identified. These liquids are characterized by hav- 
ing strong correlation in the NVT ensemble between the 
equilibrium fluctuations of the instantaneous potential 
energy U and the virial W . Recall that the instanta- 
neous energy E and pressure p can be written as a sum 
of a kinetic part and a configurational part: E = K + U 
and pV — NksT + W, respectively. The correlation be- 
tween U and W is quantified via the linear correlation 
coefficient R defined as 



R 



{AWAU) 
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correlation coefficient R = 1, since AW — {n/3)AU, and 
only IPL systems are perfectly correlating. In the study 
of strongly correlating liquids it was discovered that they 
obey Rosenfeld's excess entropy scaling, isochronal su- 
perposition, as well as density scaling. These types of 
scalings can be explained in the framework of so-called 
isomorphs (definition follows later). 



Model systems that have been identifie d "'^^1^^'^^" — 
to belong to this class of liquids include the stan- 
dard single-component Lennard- Jones liquid (SCLJ), the 
Kob- Andersen binary LJ mixture^ii^ {KABLJ), the 
asymmetric dumbbell modeP°, the Lewis- Wahnstrom o- 
terphenyl mode l"Si . ? . 'i (OTP), and several others. Strong 
WU correlation has been experimentally verified for a 
molecular van der Waals liquid^^ and for supercritical 
argonii. The class of strongly correlating liquids is be- 
lieved to include most or all van der Waals and metal- 
lic liquids, whereas covalently, hydrogen-bonding or ionic 
liquids are not strongly correlating^^. The latter refiects 
the fact that competing interactions tend to destroy the 
strong correlation. 



The class of strongly correlating liquids is defined by 
R > 0.90^^. An inverse power-law (IPL) system r~" has 



An example of strong W^C/-correlation is given in Fig. 
[T] for the asymmetric dumbbell model^*' (details of this 
model are provided in Sec. IIIip . Figure [IJa) shows the 
time evolution of the equilibrium fluctuations of U and 
W normalized to zero mean and unity standard devia- 
tion. Fig. [TJb) shows a scatter plot of the corresponding 
instantaneous values of U and W. U and W are clearly 
strongly correlated in their equilibrium fluctuations. 
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FIG. 1. Two difTerent ways of visualizing the strong 
virial/potential energy correlation for the asymmetric dumb- 
bell model at p = 1.863 and T = 0.465 (see Sec. [TTT]for details 
of the model and the units used), (a) The time evolution of 
U (black) and W (red) per particle normalized to zero mean 
and unity standard deviation, (b) A scatter plot of the in- 
stantaneous values of W and U per particle. The correlation 
coefficient R is 0.96. 

References [H and [l^ identified the cause of strong 
W^f/-correlation in the SCLJ Hquid. The LJ pair po- 
tential can be weh approximated from about r — 0.95cr 
to r = 1.5cr (see Pedersen et al3^) by a sum of an IPL, 
a linear term, and a constant via the so-called " extended 
IPL potential"i2,: VLjir) ~ Ar'"- + B + Cr. At mod- 
erate pressures this covers the entire first peak of the 
radial distribution function, i.e., the first coordination 
shell. The constraint of constant volume in the NVT 
ensemble has the effect that when one nearest neighbor 
distance increases another one decreases; upon summa- 
tion the contribution from the linear term to U and W 
is almost constant. The latter observation has the con- 
sequence that some of the scaling properties of pure IPL 
systems are inherited in the LJ system in the form of 
isomorphs. 

Reference [l^ introduced a new concept referring to a 
strongly correlating liquid's phase diagram, namely iso- 
morphic curves or more briefly: isomorphs. Two state 
points with density and temperature (pi, Ti) and {p2, 
T2) are defined to be isomorphio^i if the following holds: 
Whenever a configuration of state point (1) and one 



of state point (2) have the same reduced coordinates 
(p^^'^i'l^-' = py^r)^'^ for all particles i), these two con- 
figurations have proportional Boltzmann factors, i.e.. 

Here C12 is a constant that depends only on the state 
points (1) and (2). An isomorph is defined as a continu- 
ous curve of state points that are all pairwise isomorphic. 
In other words, Eq. ^ defines an equivalence relation 
with the equivalence classes being the isomorphs. Only 
IPL systems have exact isomorphs; these are character- 
ized by having /T = const where 7 = n/3. Reference 
[T5I argued and demonstrated by simulations that strongly 
correlating liquids have isomorphs to a good approxima- 
tion. 

From the defining property of an isomorph [Eq. ([2])] it 
follows that the structure in reduced units (f; = p^/'^Yi) 
is invariant along an isomorph, since the proportional- 
ity constant C12 disappears when normahzing the con- 
figurational canonical probabilities^^. Thus the reduced 
unit radial distribution function and the excess entropy 
Sex = S" — Sid are isomorph invariants, where Sid is the 
ideal gas contribution to the entropy at the same temper- 
ature and density. Isomorph invariance is, however, not 
limited to static quantites, also the mean-square displace- 
ment, time auto-correlation functions, and higher-order 
correlation functions are invariant in reduced units along 
an isomorph. The reader is referred to Ref. [l^ for a de- 
tailed description of isomorph invariants, and the proof 
that a liquid is strongly correlating if and only if it has 
good isomorphs. A brief overview of strongly correlating 
liquids and their isomorphs can be found in Pedersen et 
a/26. 

Reference jT^ studied isomorphs of atomic single- and 
multicomponent LJ liquids with generalized exponents 
m and n. It was found that for given exponents (m, n) 
all isomorphs have the same shape in the WU -phase dia- 
gram, i.e., a so-called master isomorph exists from which 
all isomorphs can be generated via a simple scaling of the 
VFJ7-coordinates. For instance, the shape of isomorphs 
in the W^[/-phase diagram of the SCLJ liquid and the 
KABLJ liquid are the same. 

References [1211161 focused on understanding strong 
M^?7-correlation and its implication for atomic systems. 
Schr0der et al^ in 2008 studied two rigid molecular liq- 
uids that are strongly correlating: the asymmetric dumb- 
bell model and the Lewis- Wahnstrom OTP model (see 
Sec. imp . At that time the isomorph concept had not yet 
been developed and state points with the same p^ /T , as 
inspired from the IPL system, were tested for collapse of, 
for instance, the reduced unit radial distribution function 
(note that in Refs. [l^, [H, and [13 7 is defined slightly 
different from subsequent papers). The dynamics in re- 
duced units was also found to be a function of p"^ jT ^ to 
a good approximation, as is the case for IPL systemsii. 
Chopra et al^ found that the Sex can be written (ap- 
proximately) as a function of p'^ /T for rigid symmetric 
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LJ dumbbells with different bond lengths. They also 
found that the reduced diffusion constant and relaxation 
time are functions of Sex- These results suggest that 
the isomorph concept is relevant also for rigid molecular 
systems. In this paper we expand on earlier results by 
studying in detail the same systems as Schr0der et al?^ . 

The isomorph definition Eq. ^ must be modified for 
rigid molecules, since the bond lengths are fixed and can- 
not follow the overall scaling. A simple modification of 
Eq. ([2]), which is consistent with the atomic definition, 
is to define the mapping amongst configurations in terms 
of the molecular center of masses, instead of the atomic 
positions. We thus define two state points in the phase 
diagram of a liquid composed of rigid molecules to be 
isomorphic if the following holds: Whenever two configu- 
rations of the state points have identical reduced center- 
of-mass coordinates for all molecules, 



1/3 (1) _ 1/3 (2) 
Pi ^CM,i — P2 ^CM,i^ 



as well as identical Eulerian angle; 
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(4) 



these two configurations have proportional Boltzmann 
factors, i.e., [where R = (rcj\/,i, di, ^i, i/'i, rcAf.w, 
On, i'N, V'A')] 



-C/(R(i')/fcBTi ^ (^^^g-C/(R(2))/fcsT2^ 



(5) 



Again, C12 is a constant that depends only on the state 
points (1) and (2). An isomorph is defined as a set of 
state points that are pairwise isomorphic. It should be 
noted that, in contrast to what is the case for atomic sys- 
tems, since the bonds do not follow the overall scaling of 
the system this definition does not imply the existence of 
exact isomorphs for rigid molecules with intermolecular 
IPL interactions,. 

Taking the logarithm of Eq. ([5]) implies 



i7(R^^0 - 7^2/^1 ■ C/(R^^O + kBT2 In C12. (6) 

Equation ([5]) provides a convenient way of testing to 
which extent Eq. ([5]) is obeyed for a given system. A 
simulation is performed at one state point (1) and the ob- 
tained configurations are scaled to a different density p2, 
where the potential energy is evaluated. The respective 
potential energies of the two state points are then plotted 
against each other. In the resulting plot a near straight- 
line indicates that there exists an isomorphic state point 
with density p2- The temperature T2 of the isomorphic 
state point can be found from the slope of a linear re- 
gression fit. This procedure is termed the "direct iso- 
morph check"—. If this test is performed for an atomic 
IPL system a correlation coefficient of i? = 1 is obtained, 
consistent with these systems having exact isomorphs. 



As an example we perform a direct isomorph check for 
the asymmetric dumbbell model in Fig. [21 A correlation 
coefficient oi R = 0.97 is observed for a 15% density 
increase. Calculating the temperature of the isomorphic 
state point from the linear regression slope the result 
differs only 1% from the prediction by requiring constant 
excess entropy (see Sec. HV]). 
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FIG. 2. The "direct isomorph check""iS for the asymmetric 
dumbbell model. During a simulation at state point (pi, Ti) 
= (1.737, 0.309) the center of mass of each dumbbell is scaled 
to density p2 = 1.998, keeping the Eulerian angles fixed. The 
potential energy is then evaluated from the scaled configura- 
tions and plotted against the potential energy of the unsealed 
configurations. The temperature r2 (slope) of the isomorphic 
state point at density p2 is calculated by multiplying the lin- 
ear regression slope with Ti [Eq. ([S])]. 

In the present paper we show that liquids composed of 
simple rigid molecules have good isomorphs in their phase 
diagram as defined in Eqs. ([3]) - ([S]). Section HIl derives 
several isomorph invariants in molecular systems com- 
posed of rigid molecules. Section IIIII describes the sim- 
ulation setup and the investigated model systems. Sec- 
tion IIVI investigates the existence of isomorphs for the 
asymmetric dumbbell and the Lewis- Wahnstrom OTP 
models^i^i. Section |V] investigates the existence of a 
master isomorphic for these model systems. Section IVII 
summarizes the results and presents an outlook. 



II. ISOMORPH INVARIANTS IN LIQUIDS 
COMPOSED OF RIGID MOLECULES 

From the single assumption of curves of isomorphic 
state points in an atomic liquid's phase diagram, Ref. [l^ 
derived several invariants along an isomorph. Since we 
have extended this definition in Eqs. ([3]) - (O to molec- 
ular systems composed of rigid molecules, it is natural 
to wonder which of these invariants can be extended to 
molecular systems. The molecular isomorph concept is 
different from the atomic case in that there is no " ideal" 
reference system (the IPL system). The following sec- 
tions and simulations, however, show that isomorphs can 
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nevertheless be a useful tool for understanding such liq- 
uids. 

In the following we derive several invariants from exact 
isomorphs. We start by noting that the generalization of 
isomorphs to molecular systems define a bijective map 
amongst configurations of state points (1) and (2). The 
NVT configurational probability density for a system of 
N rigid molecules is given by^^ (where dR = dr^j^jdr^ 
with T = {9, (f), ip) and dr — sin ddd(j)dip) 



P(R) 



-U(-R.)/kBT 



Jg-C/(R)/fcBT^R- 



(7) 



In combination with Eq. ([5|) it follows that all mapped 
configurations of state points (1) and (2) have identical 
Boltzmann probabilities, i.e., 



(2) 



(8) 



For convenience we introduce two configurational distri- 
bution functionsi^iSS 



P(R) = {VnYP{Il), 



(9) 

, (10) 



where is the integral over the Eulerian angles for 
one molecule {fl = Stt^ for a non-linear molecule). 
P has been introduced to make P dimensionless. 
P{r^n[,T^)dr^i^dT^ is the probability to observe the 
system represented by a point in the volume-element 
dr^j^^jdr^ located at {fcMi'''^}- ^ invariant along 
an isomorph and is related to P via P(fp^^,T^) = 
{Nn)-'^P{TL) = p-^P(R). We note that the excess 
entropy Sex = —{dFex/dT)iqy, where F^x is the excess 
free energy, can be written asi^ 



Sex = -fcs j (ra)-^P(R) lnP(R)dR, (11) 
= -ks I PhiP dr^^dr^ - kgNlniNn). (12) 



From the above observations we can now derive a num- 
ber of isomorph invariants in liquids composed of rigid 
molecules. 

1. The molecular center- of-mass structure in reduced 
units. For a given configuration of the molecular 
center-of-mass structure in reduced units, all ori- 
entations of the molecules of state points (1) and 
(2) by Eq. ([8|) have identical probabilities. The 
reduced center-of-mass structure is thus invariant 
along an isomorph. 



2. Any normalized distribution function describing the 
(relative) orientations of molecules with respect 
to their center-of-mass. All orientations of the 
molecules with respect to a given molecular center- 
of-mass configuration of state point (1) are mapped 
to configurations of state point (2) with identi- 
cal probabilities. It thus follows that the normal- 
ized distribution function is invariant along an iso- 
morph. 



3. The isochoric heat capacity Cv ■ The excess heat 
capacity in the NVT ensemble is given by Cy = 
{{AUy)/kBT^. Defining X = U/ksT we may 
write C^^ = fcB((AX)2). By Eqs. (0) and © it 
follows that Cy^ is invariant along an isomorph, 
since the constant kBT2 In C12 disappears when 
subtracting the mean. The ideal gas contribution 
to Cv is independent of state point {Cy = 6NkB/2 
for non- linear molecule). 

4. The translational two-body entropy^^i^^^ St/N = 
-pkB/2 J[gcM{r) In gcMir)-gcMir)-\-l]dr, where 
gcMif) is the radial distribution function for the 
center of mass of the molecules. The density depen- 
dence disappears when switching to reduced units 
and by Statement [T] the molecular center-of-mass 
structure in reduced units is invariant along an iso- 
morph, and thus also the radial distribution func- 
tion (in reduced units). 

5. The orientational two-body entropy^^/^^^ So/N 
- -pfcs/(2f]2) / gcM{r)9{co^\r)\ng{cu^\r)doj^dr, 
where oj^ denotes a set of angles used to describe 
the relative orientation of two molecules, and 
g(a;^|r) is the conditional distribution function for 
the relative orientation of two molecules separated 
by a distance r. Applying reduced units this 
invariant follows from Statements [TJ and [5] 

6. All N-body entropy terms^^^. The excess entropy 
can be expanded as Sex = Si^i ^i- The two-body 
expression S2 = St -\- So is given above, where as 
the higher-order terms are more involved. 

7. The excess entropy Sex ■ The excess entropy is given 
by Eq. (|12p . and since P is invariant along iso- 
morph, so is the excess entropy. The latter also 
follows from Statement [6l, since each term is in- 
variant. 

8. The molecular center-of-mass NVE and Nose- 
Hoover NVT dynamics in reduced units. The re- 
duced dynamics of the atomic positions on account 
of the constraints is not invariant along an iso- 
morph. Considering instead the molecular center- 
of-mass motion the constraint force disappears and 
these equations of motion are invariant along an 
isomorph in reduced units. The proof is given in 
Appendix |B] (a brief summary of constrained dy- 
namics is given in Appendix IXj) . 
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9. Any average molecular center-of-mass dynamic 
quantity in reduced units. This follows immedi- 
ately from Statement m, since the molecular center- 
of-mass equations of motion in reduced units are 
invariant along an isomorph. In particular, this 
would include a reduced relaxation time fa- 

As detailed above it is necessary to consider the (reduced) 
center of mass motion and the motion relative to the cen- 
ter of mass separately. Nevertheless, during the investi- 
gation of isomorphs we will also consider the reduced 
atomic quantities to examine their " invariance" . 

An additional consequence of isomorphs is that, since 
by Eq. two isomorphic state points have identical 
canonical probabilities, an instantaneous change of tem- 
perature and density from an equilibrated state point to 
another isomorphic state point does not lead to any re- 
laxation. This is called an isomorphic jumpi^. 



III. SIMULATION DETAILS 



The potential energy has no contribution from the fixed 
bonds. A force smoothing procedure22^ was applied from 
rs = 2.45(Tq,^ to Tc — 2.50(Tq^. 

The bond lengths were held fixed using the Time 
Symmetrical Central Difference (TSCD) algorith m^**'^^ , 
which is a central difference time-discretization of 
the constrained equations of motion preserving time- 
reversibility. The simulations were performed in the NVT 
ensemble applying the Nose-Hoover {NH) algorithmic- 
using R UMD^, a molecular dynamics package optimized 
for state of the art GPJJ-computing. 

The NVT simulations were performed without adjust- 
ing the time-constant of the NH algorithm (see Appendix 
|B|) . This choice is not expected to influence the results 
over the observed density and temperature range, since 
the dynamics is not particular sensitive to the absolute 
value of the NH time-constant^. 

Appendix |X] gives a brief summary of constrained dy- 
namics and the effect on the virial (see also Refs. ISJ, 
[ill , and [42! ) . The specific details of the investigated two 
models follow below. 



We studied two model systems of rigid molecules (Fig. 
[3]): the asymmetric dumbbell model {N = 500) and the 
Lewis- Wahnstrom OTP model {N = 320). 
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The asymmetric dumbbell Lewis-Wahnstrom OTP 



FIG. 3. A sketch of the two model systems studied: The 
asymmetric dumbbell and Lewis-Wahnstrom OTP models. 
The asymmetric dumbbell is a simplistic model of toluene 
with the methyl side-group tightly bonded to the benzene 
molecule. The Lewis-Wahnstrom OTP model is an isosceles 
triangle with an angle of 75°, different from the 60° of the 
real 1,2-diphenylbenzene molecule^^. 

For both models the potential energy U and the virial 
W are given by 



A. The asymmetric dumbbell 

The asymmetric dumbbell model consists of a big (A) 
and small (B) LJ particle, rigidly bonded with bond dis- 
tance of rij = 0.584 (here and henceforth units are given 
in LJ units referring to the A particle, oaa = 1, ^aa = 
1, and tua = !)• The parameters were chosen to roughly 
mimic Toluene^S. The asymmetric dumbbell model can 
be cooled to a highly viscous state without crystallizing, 
making it feasible to study slow dynamics. The asym- 
metric dumbbell model has cfab — 0.894, (Jbb = 0.788, 
tAB = 0.342, and €bb = 0.117 with 0.195. 



B. Lewis-Wahnstrom OTP 

The Lewis-Wahnstrom OTP model^^ii^ consists of 
three identical LJ particles rigidly bonded together in 
an isosceles triangle with sides of = 1.000 and top- 
angle of 75°, i.e., different from the 60° of the real 1,2- 
diphenylbenzene molecule^. All parameters (including 
the masses) are unity for the OTP model. 



W ^Wlj + Wcon- 



(13) 
(14) 



The first term in the virial is the LJ virial Wlj, the 
second term is the contribution to the virial due to the 
constraints (fixed bond lengths), Wcon- Ulj is a sum 
over intermolecular pair interactions given by the (12, 
6)-LJ potential 



(15) 



IV. NUMERICAL STUDY OF ISOMORPHS 
FOR THE TWO MODEL SYSTEMS 

In order to investigate whether the two model systems 
have good isomorphs, we first describe how to generate 
an isomorph in a simulation. The excess entropy Sex is 
invariant along an isomorph, and a method for generating 
an isomorph is to generate a curve of constant Sex (see 
Sec. im and also Refs. [l^ and flGh. A curve of constant 
excess entropy can be found by using the exact NVT 
ensemble relatiorti^ 



(A[/AVF) 



i9 In p / 



(16) 



In simulations an isomorph is generated as follows: 1) 
The left-hand side is calculated from the fluctuations at 
a given state point; 2) A new state point is identified 
by a dicretization of Eq. ([T6| by changing the density 
by 1%, where the new temperature is calculated from 
AlnT = 7Alnp; 3) The procedure is repeated and in 
this way an isomorph is generated in the phase diagram. 
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A. Isomorphs of the asymmetric dumbbell model 



This section investigates the asymmetric dumbbell 
model. Isomorphs were mapped out following the pro- 
cedure described above. Figure 2] shows the AA radial 
distribution functions along an isomorph with 19% den- 
sity increase before (a) and after (b) scaling the distance 
r into reduced units 
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Also shown for reference in Fig. |4jc) is the AA radial 
distribution functions along an isotherm with 12% 
density increase. Although the reduced structure of the 
atomic positions is not predicted to be invariant along 
an isomorph, Fig. |4] shows that it nevertheless is a 
resonable approximation. The reduced structure of the 
atomic positions is less invariant along the isotherm. 



FIG. 4. Radial distribution functions for the asymmetric 
dumbbell model, (a) AA pair-correlation function along an 
isomorph with 19% density increase before scaling the dis- 
tance r into reduced units f = p^^^r. (b) AA pair-correlation 
function along the same isomorph after scaling the distance r 
into reduced units, (c) AA pair-correlation function along an 
isotherm with 12% density increase function after scaling of 
the distance. 



Figure [5] considers the AB radial distribution func- 
tions, where the constrained bond distance shows up as 
a sharp peak. The analogous conclusion as with the AA 
distribution functions is reached, and likewise for the BB 
distribution functions fnot shown). 
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FIG. 5. Radial distribution functions for the asymmetric 
dumbbell model, (a) AB pair-correlation function along the 
isomorph of Fig. U before scaling the distance r into reduced 



units f = p 



1/3^ 



(b) AB pair-correlation function along the 



same isomorph after scaling the distance r into reduced units, 
(c) AB pair-correlation function along the isotherm of Fig. |4] 
after scaling of the distance. 



Next, we consider in Fig. [6] the molecular center-of- 
mass radial distribution functions along the isomorph 
and isotherm of Figs. HUH This quantity is predicted 
to be invariant along an isomorph (see Sec. |lTl. The 
molecular center-of-mass structure is to a good approx- 
imation invariant in reduced units along the isomorph, 
while this is less so along the isotherm as can be seen 
from the first peak. 
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FIG. 6. Molecular center-of-mass radial distribution func- 
tions for the asymmetric dumbbell model, (a) Pair-correlation 
function along the isomorph of Figs. |4][5] before scaling the 
distance r into reduced units f = p^^^r. (b) Pair-correlation 
function along the same isomorph after scaling the distance 
r into reduced units, (c) Pair-correlation function along the 
isotherm of Figs. |4][5] after scaling of the distance. 



We consider in Fig. [7] the dynamics in terms of the re- 
duced A-particle incoherent intermediate scattering func- 
tion. The reduced dynamics of the atoms is not predicted 
to be invariant along an isomorph (see Appendix[B]) , how- 
ever, the figure shows that it is a good approximation. 
The same conclusion is reached for the B-particle (not 
sho-w" 
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FIG. 7. The reduced yl-particle incoherent intermediate scat- 
tering function for the asymmetric dumbbell model keeping 
the reduced wavevector q constant, (a) Along the isotherm of 
Figs. |3]iniwith 12% density increase, (b) Along the isomorph 
of Figs. 1 4161 with 19% density increase. 



FIG. 8. The reduced molecular center-of-mass incoherent in- 
termediate scattering function for the asymmetric dumbbell 
model keeping the reduced wavevector q constant, (a) Along 
the isotherm of Figs. I4I7I with 12% density increase, (b) 
Along the isomorph of Figs. 141 71 with 19% density increase. 



In Figure [5] we consider the reduced molecular center- 
of-mass self part of the intermediate scattering function. 
This quantity is predicted to be invariant along an 
isomorph (see Appendix [B]), and Fig. |8] clearly shows 
this. The dynamics is not invariant along the isotherm. 



We show the variation of 7, calculated from the NVT 
fluctuations via Eq. (fTB)) . in Fig. [9] along an isochore 
and along the isomorph of Figs. |3][5]in two different ver- 
sions. The crosses show 7 calculated from the total virial 
W while the asterisks show 7 calculated after subtract- 
ing the constraint contribution to virial, i.e., replacing 
W with Wlj — W — WcoN- The insets show the cor- 
responding correlation coefficients R. Reference [l^ pre- 
dicts that 7 should be a function of density only 7 = 
7(/9). This is seen in Fig. [9] to be a good approximation 
for both versions of 7, where the crosses are the 7 used 
to keep the excess entropy constant. 
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As mentioned in the introduction, density scalingir"— 
is the empirical observation that the relaxation time Tq 
for many viscous liquids can be written as some function 
'''a = f{p'^''"^^'/T) where ^scaie in experiments is a fitting 
exponent. In performing these fits reduced units are 
often not used; however, the importance of using re- 
duced units in experimental data has only recently been 
pointed out*'^. If we assume that 7 is constant along an 
isomorph, Eq. implies that p^' /T = const describes 
the isomorph. In this case density scaling will thus hold 
to a good approximation since the reduced relaxation 
time is an isomorph invariant^''; for the dumbbell system 
7 changes only moderately along an isomorph and thus 
density scaling is a good approximation for this system. 
That 7 for systems with isomorphs can be identified with 
the density scaling exponent ^scaie has very recently 
been verified expementially for a silicone oil^^. 
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FIG. 9. The variation of 7 [Eq. (|16|l ] and the correlation 
coefficient 7? [Eq. ([1}] for the asymmetric dumbbell model in 
two different versions along an isochore (red, p — 0.932) and 
along the isomorph (black) of Figs. I4I8I The crosses show 
7 calculated from the total virial W while the asterisks show 
7 calculated after subtracting the constraint contribution to 
virial, i.e., Wlj = W — Wcon- The corresponding R's are 
shown in the insets. 7 is predicted in Ref. [T^ to be a function 
only of density which is seen to apply to a good approximation 
for both versions. 
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Starting from an equilibrated sample at some state 
point, changing either temperature or density alters the 
equilibrium Boltzmann distribution of states. Two iso- 
morphic state points have identical Boltzmann canonical 
probabilities (Eq. ([5|)). A sudden change of state from 
one state point to another isomorphic state point should 
thus not lead to any relaxation. This is called an iso- 
morphic jump, and the prediction of no relaxation was 
shown in Ref. JJ, to work well for the KABLJ liquid. 



A similar numerical experiment is carried out for the 
asymmetric dumbbell model in Fig. 1101 Considering 
three equilibrated, isochoric state points (1), (2), and (3), 
the density and temperature are instantaneously changed 
to a state point (4). State point (4) is isomorphic to state 
point (2). The isomorph prediction is that jumps from 
state points (1) and (3) show relaxation, but not jumps 
from state point (2). This is indeed the case (Fig. [rUT a)). 
State point (1) ages from below since the aging scheme 
(1) (4) can be described as first an instantaneous iso- 
morphic jump to the correct density, but a lower temper- 
ature, and subsequently relaxation from this state point 
along the isochore of state point (4). 



FIG. 10. Four state points (1), (2), (3), and (4) corresponding 
to, respectively, (p, T) = (0.932, 0.400), (0.932, 0.465), (0.932, 
2.000), and (0.851, 0.274) are given where the first three state 
points are on the same isochore. State points (2) and (4) are 
isomorphic, while (1) and (3) are not isomorphic to (4). After 
equilibrating at state points (1), (2), and (3), respectively, the 
temperature and density are instantaneously changed to that 
of state point (4) via a scaling of the coordinates keeping 
bond distances and Eulerian angles of the molecules fixed. 
An average has been performed over 100 samples, (a) The 
relaxational behaviour of all state points quantified by the 
potential energy U . The jump (2) — >■ (4) shows no relaxation 
while the other state points do. (b) Close up of the potential 
energy of state point (2) before and after the jump, where the 
jump takes place at f « 60. 



We finally consider the excess heat capacity per parti- 
cle Cy /N in Fig. [TT] along the isomorph and isotherm 
of Figs. I31151 The excess heat capacity increases less than 
2% along the isomorph, while the 12% density increase 
on the isotherm results in a 25% increase in the excess 
heat capacity. This is consistent with the prediction in 
Sec. |TT]that /N is an isomorph invariant. 
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FIG. 11. The excess heat capacity per particle /N for the 
asymmetric dumbbeU model as a function of density along the 
isomorph (black) and isotherm (red, T — 0.465) of Figs. I4I8I 
The density increase is 19% and 12%, respectively. The excess 
heat capacity increases less than 2% along the isomorph, while 
the isotherm shows a 25% increase. For the isotherm the 
dynamics becomes very slow for densities higher than p — 
0.950 and the system becomes difficult to equilibrate properly. 
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FIG. 12. The correlation coefficient 7? as a function of the 
bond length in the asymmetric dumbbell model at (p, T) 
= (0.932, 0.465). The investigated model has bond length 
0.584 with a correlation coefficient R « 0.97; however, as 
the bond length increases, the correlation coefficient decreases 
to i? ~ 0.65 at unity bond length. The insets shows the 
corresponding values of 7 as defined in Eq. p6|) . As the 
bond length increase the system becomes very viscous and 
the statistics is poor at high bond lengths. 



B. Isomorphs of a symmetric IPL dumbbell model 



The previous figures show that isomorphs exists to a 
good approximation for the asymmetric dumbbell model. 
An important question is whether the specific molecu- 
lar geometry determines whether or not a particular LJ 
model has good isomorphs. In Fig. [12] the correlation 
coefficient R is given as a function of the bond length. 
The correlation coefficient decreases to i? « 0.65 at unity 
bond length, and thus one might be tempted to conclude 
that LJ models with large bonds lengths in general do 
not have good isomorphs. In Sec. IIV Cl we investigate the 
Lewis- Wahnstrom OTP model which have unity bond 
lengths and show that this model actually has good iso- 
morphs. A theory connecting the variation of R to the 
molecular geometry and/or bond lengths remains to be 
developed. 



In this section we briefly consider the isomorphs of a 
symmetric inverse power-law (IPL) dumbbell with ex- 
ponent n = 18 and bond length 0.584. In Figs. [T4lja) 
and (b) we show the particle radial distribution functions 
along an isomorph before and after scaling the distance 
r into reduced units. Also shown is the reduced particle 
incoherent intermediate scattering function in Fig. 114( c). 
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FIG. 13. Structure and dynamics along an isomorph with 
19% density increase for the symmetric IPL dumbbell model 
(n = 18). (a) Particle pair-correlation function before scaling 
the distance r into reduced units, (b) Particle pair-correlation 
function after scaling the distance r into reduced units, (c) 
The reduced particle incoherent intermediate scattering func- 
tion at constant reduced wavevector q. 
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The corresponding molecular center-of-mass quantities 
are shown in Fig. [TH Interestingly, the atomic dynamics 
appears more invariant than for the reduced molecular 
dynamics. The latter is predicted to be invariant along 
an isomorph while the former is not; however, we have 
not tried to quantify this observation any further. 

Atomic systems with IPL interactions have exact 
isomorphs. This reflects the scale invariance of the IPL 
potential, i.e., that it preserves its shape under a scaling 
of the argument. Since molecules by their fixed geometry 
define a length scale in the system, isomorphs will always 
be approximate. However, the previous figures show 
that rigid molecules with IPL intermolecular interactions 
can also have good isomorphs. 



FIG. 14. Structure and dynamics along the isomorph of Fig. 
1131 for the symmetric IPL dumbbell model (n = 18). (a) 
Molecular center-of-mass pair-correlation function before scal- 
ing the distance r into reduced units r = p^^^r. (b) Molecu- 
lar center-of-mass pair-correlation function after scaling the 
distance r into reduced units, (c) The reduced molecular 
center-of-mass incoherent intermediate scattering function at 
constant reduced wavevector q. 

In Fig. [15] we consider the variation of 7 and R along 
the isomorph, which shows that R decreases slightly with 
increasing temperature and density. The variation of 7 
along the isomorph is less than for the asymmetric dumb- 
bell, and is to a good approximation constant. As for 
the asymmetric dumbbell model (Fig. [5]) the effect of 
the constraints is to increase 7 and decrease R (these are 
respectively 6 and 1 for the IPL potential used). 
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FIG. 15. The variation of 7 and R (insets) along the iso- 
morph of Figs. I13I14I with 19% density increase for a sym- 
metric dumbbell model (n = 18). 7 increases slightly along 
the isomorph. Excluding the constraints in the virial the cor- 
relation coefficient and 7 are respectively R = 1 and 7 = 6. 
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We proceed to investigate the Lewis— Wahnstrom OTP 
mode l '— 1^ '. Figure [TBI shows the particle radial distribu- 
tion functions along an isomorph with 21% density in- 
crease before and after scaling the distance r into reduced 
units. We treat the particles as identical in the quantities 
probed in simulations (i.e., the radial distribution func- 
tion, etc.) even though the OTP model is an isosceles 
triangle. Also shown for reference is an isotherm with 
11% density increase in Fig. IWc). 



FIG. 16. Particle radial distribution functions for the OTP 
model, (a) Along an isomorph with 21% density increase, 
shown prior to scaling the distance r into reduced unit f — 
p^^'^r. (b) Along the same isomorph after scaling the distance 
r into reduced unit, (c) Along an isotherm with 11% den- 
sity increase. At the highest density probed the OTP model 
crystallizea^^ (magenta). 
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Figure [T7] shows the corresponding reduced molecular 
center-of-mass radial distribution functions. The reduced 
molecular center-of-mass structure seems less invariant 
along the isomorph than for the asymmetric dumbbell 
(Fig. [5]), consistent with the OTP model being less 
strongly correlating (i? « 0.90). However, comparing 
with the isotherm in Fig. [TTT c) the OTP model crys- 
tallizes at the highest density probed"^^ , even though the 
density increase is just 11% compared with the 21% den- 
sity increase along the isomorph. Comparing now with 
the particle quantities of Fig. [16] the latter seems more 
invariant, even though the reduced molecular center-of- 
mass structure is predicted to be invariant. We currently 
have no explanation for this observation. 
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FIG. 18. The reduced particle incoherent intermediate scat- 
tering functions for the OTP model keeping the reduced 
wavevector q constant, (a) Along the isotherm of Figs. 11611171 
with 11% density increase, (b) Along the isomorph of Figs. 
I16I17I with 21% density increase and almost a factor of four 
increase in temperature. The dynamics is roughly invariant 
along the isomorph, but not along the isotherm. 



For the molecular quantities, the dynamics is roughly 
invariant along the isomorph but not on the isotherm, 
even though the density increase is 21% for the isomorph 
and only 11% for the isotherm. In contrast to the re- 
duced molecular center-of-mass structure; the molecular 
dynamics does not seem less invariant than the particle 
dynamics, consistent with the prediction of Appendix IB] 



FIG. 17. Molecular center-of-mass radial distribution func- 
tions for the OTP model, (a) Along the isomorph of Fig. [16] 
with 21% density increase, shown prior to scaling the distance 
r into reduced units f = p^^^r. (b) Along the same isomorph 
after scaling the distance r into reduced units, (c) Along the 
isotherm of Fig. [T6]with 11% density increase. At the highest 
density probed the OTP model here crystallizes (magenta). 



Figure [TH] shows the reduced particle incoherent in- 
termediate scattering functions along the isotherm and 
isomorph of Figs. [TBlfTTl while Fig. [12] shows the re- 
duced molecular center-of-mass incoherent intermediate 
scattering functions. 
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FIG. 19. The reduced molecular center-of-mass incoherent 
intermediate scattering functions for the OTP model keeping 
the reduced wavevector q constant, (a) Along the isotherm 
of Figs. 11611181 with 11% density increase, (b) Along the 
isomorph of Figs. 11611181 with 21% density increase and al- 
most a factor of four increase in temperature. The dynamics 
is roughly invariant along the isomorph, but not along the 
isotherm. 



We consider in Fig. [5D]the variation of 7 as defined by 
Eq. ([TC)) . The large variation in 7 indicates that density 
scaling, with a fixed exponent, is of more approximate 
nature for the OTP modei^S than for the asymmetric 
dumbbell, even though isomorphs exists to a good 
approximation. The isomorphs of the OTP model are, 
however, more approximative than for the asymmetric 
dumbbell, which is consistent with OTP model being 
less strongly correlating. 
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FIG. 20. The variation of 7 (Eq. (|16|) ) and the correlation co- 
efficient R (inset, Eq. ((T|)) for the OTP model in two different 
versions along an isochore (red, p = 0.329) and the isomorph 
(black) of Figs. I16I19I The crosses show 7 calculated from 
the total virial W , the asterisks show 7 calculated after sub- 
tracting the constraint contribution to virial, i.e., replacing 
W by Wlj = W — WcoN- 7 is predicted in Ref. 15 to be a 
function only of density as is seen to be the case for both ver- 
sions, although the variation is larger than for the asymmetric 
dumbbell. 

Next we consider isomorph jumps for the OTP model. 



The setup is analogous to that of the asymmetric 
dumbbell model described in Sec. IIV Al It is seen from 
Fig. [5T]that an isomorph jump shows no relaxation. 
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FIG. 21. Four state points (1), (2), (3), and (4) corresponding 
to, respectively, (p, T) = (0.329, 0.650), (0.329, 0.700), (0.329, 
1.000), and (0.303, 0.383) are given where the first three state 
points are on the same isochore. State points (2) and (4) are 
isomorphic while (1) and (3) are not isomorphic to (4). After 
equilibrating at state points (1), (2), and (3), respectively, 
the temperature and density are instantaneously changed to 
that of state point (4) via a scaling of the coordinates keeping 
bond distances and Eulerian angles of the molecules fixed. 
An average has been performed over 100 samples, (a) The 
relaxational behaviour of all state points quantified by the 
potential energy U. The jump (2) — >■ (4) shows no relaxation 
while the other jumps do. (b) Close up of the potential energy 
of state point (2) before and after the jump, where the jump 
takes place at f ~ 60. 



We close the investigation of the OTP model by con- 
sidering in Fig. [22] the excess heat capacity per particle 
/N. This quantity increases 7% over the 21% den- 
sity increase along the isomorph, while the 11% density 
increase on the isotherm results in a 34% increase be- 
fore crystallizing. These results are consistent with the 
prediction that Cy^ is an isomorph invariant (see Sec. 
HI]) , although less so than for the asymmetric dumbbell 
model. 
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FIG. 22. The excess heat capacity per particle C'^ /N for 
the OTP model as a function of density along the isomorph 
(black) and isotherm (red) of Figs. 11611191 The density in- 
crease is 21% and 11%, respectively. At high densities for the 
isotherm the OTP model crystallizes. The excess heat capac- 
ity is to a good approximation invariant along the isomorph, 
while this is not the case for the isotherm. 
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FIG. 23. (a) Three different isomorphs for the asymmet- 
ric dumbbell model in two different versions with 19%, 21% 
and 22% density increase, respectively (black, magenta and 
green). The crosses give the total virial W, the asterisks give 
Wlj = W — WcoN- fa is reduced relax;ation time extracted 
from the self part of the intermediate scattering function, (b) 
The same isomorphs as in (a) where WU and WljU are scaled 
to superpose with a factor identified by trial and error. The 
black points have unity scaling factor. 



MASTER ISOMORPHS 



The previous section detailed the existence of iso- 
morphs in the phase diagram of liquids of small rigid 
molecules. We now investigate whether the generated 
isomorphs have the same shape in the WU phase dia- 
gram, i.e., whether a so-called master isomorph existsi^. 
It is also interesting to compare the isomorphs of the 
dumbbell and the OTP models, since both systems have 
intermolecular (12, 6)-LJ interactions, but different con- 
straint contributions to the virial (one versus three con- 
strained distances per molecule). 

Figure [23l a) shows three different isomorphs in the 
WU phase diagram for the asymmetric dumbbell model, 
in two different versions: one for the total virial W 
and one for the "LJ" virial, i.e., replacing W with 
Wlj — W — WcoN- In order to investigate whether 
a master isomorph exists Fig. I^STb) shows the same 
isomorphs after scaling of the potential energy and the 
virial with the same factor (depending on the isomorph). 
The best scaling factor was identified by trial and error. 



Corresponding figures for the OTP system are given in 
Fig. [M] The figures show that for both models a master 
isomorph exists to a good approximation both with and 
without the constraint contribution to the virial. 
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FIG. 24. (a) Two different isomorphs for the OTP model 
in two different versions with 21% and 14% density increase 
(black and magenta). The crosses give the total virial W, the 
asterisks give Wlj — W — Wcon- Tc is reduced relaxation 
time extracted from the self part of the intermediate scatter- 
ing function, (b) The same isomorphs as in (a) where WU 
and WljU axe scaled to superpose with a factor identified by 
trial and error. The black points have unity scaling factor. 



As mentioned in the introduction, Ref. [Ifl derived pre- 
dictions concerning the shape of isomorphs for atomic 
systems with pair potential given by a sum of two IPLs 
(the generaUzed LJ potential). The question arises 
whether WljU follows that shape? This is studied in 
Fig- B5l a) where the WljU isomorphs for the asymmet- 
ric dumbbell and OTP models are shown scaled using 
the previously mentioned procedure. The two dashed 
curves are the isomorph prediction for an atomic (12, 6)- 
LJ system^- (where p — p/ p* and p* is the density of a 
chosen reference state point) 



u = u*^p'' + u:p\ 

LJ 



W,.j = 4t/„;p4 



2u:p\ 



(18) 
(19) 
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FIG. 25. (a) Scaled WljU isomorphs for the asymmetric 
dumbbell and OTP models. The black points have unity scal- 
ing factor. Both systems have intermolecular (12, 6)-LJ in- 
teractions, and the dashed curves are the isomorph prediction 
from Ref. [3 for an atomic system, where the LJ reference 
coefficients have respectively been calculated from the dumb- 
bell state points (p, T) = (0.932, 0.465) and (0.851, 0.274). 
(b) Scaled WU isomorphs for the systems in (a). The total 
virial does not show exact scaling between the asymmetric 
dumbbell and OTP models. 



where the reference coefficients (U^, U*) have been cal- 
culated from two different reference state points along 
"Isomorph 1" of the asymmetric dumbbell (see Ref. [la 
for details on calculating these coefficients). The only 
assumption used in Ref. [H to derive these formulas is 
the invariance of the reduced atomic structure along an 
isomorph, however, this is not predicted to be the case 
for molecular systems with isomorphs. 

It is clear that the atomic isomorph shape is not fol- 
lowed exactly. Nevertheless, there seems to exist not only 
a master isomorph in the LJ and total virial for the indi- 
vidual systems, but also for the LJ virial between these 
two different model systems. The same does not hold for 
the total virial as can been seen in Fig. [?51 b) since the 
constraint contributions are different. 



To examine the extent of "deviation" from Eqs. (fT8)) 
and we show in Fig. [51] for the asymmetric dumb- 
bell U and Wlj as a function of the reduced density 
{p* = 1). The reference coefficients can be calculated 
from a linear regression fit of the potential energy and 
the estimated coefficients can be used to plot a straight 
line in the LJ virial plot. This is performed in Fig. [^S] 
where it is clear that even though both plots follow a 
near straight-line, the coefficients are not given by Eqs. 
(IT51) and ((T5| . It is worth mentioning again that the pre- 
diction of Ref. is for an atomic system, and is as such 
not excepted to hold for rigid molecular systems. 
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FIG. 26. The potential energy and LJ virial as a function 
of for "Isomorph 1" of the asymmetric dumbbell, (a) A 
linear regression fit of the potential energy has been performed 
to calculate the reference coefficients {U^, U^). (b) These 
coefficients are then used to plot the red straight line, which 
according to the atomic prediction [Eqs. p8|l and (I19p ] should 
coincide with the black data points. The green line is a linear 
regression fit to the same data points. 



Finally we consider in Fig. [27] for the asymmetric 
dumbbell how the instantaneous fluctuations of Wcon 
correlate with Wlj and W respectively. The constraint 
contribution to the virial at this state point does not 
correlate well with the contribution to the virial com- 
ing from the LJ interactions {R = 0.31). Obviously, the 
correlation is higher when the total virial is considered 
{R = 0.61). The main contribution to the virial, for the 
asymmetric dumbbell model, comes from the LJ interac- 
tions, however, the LJ virial does not correlate well with 
the constraint virial. The latter observation may indi- 
cate a break-down of master isomorph scaling (for the 
total virial) at high pressures, however, this remains to 
be confirmed. 
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FIG. 27. (a) The correlation of the instantaneous fluctuations 
of Wlj and Wcon- The correlation coefficient 7? is 0.31. (b) 
The correlation of the instantaneous fluctuations of W and 
Wcon- The correlation coefficient R is 0.61. 



VI. SUMMARY AND OUTLOOK 

Isomorphs are curves in the phase diagram of a 
strongly correlating liquid along which a number of static 
and dynamic quantities are invariant in reduced units. 
References [151 and HQ focused on understanding iso- 
morphs in atomic systems. In this paper we general- 
ized the isomorph concept to deal with systems of rigid 
molecules [Eq. ([5])] and investigated several isomorph 
invariants for the asymmetric dumbbell model and the 
Lewis- Wahnstrom OTP model. We find that these rigid 
molecular systems also have isomorphs to a good approx- 
imation; however, the isomorphs of the OTP model were 
more approximative than those of the asymmetric dumb- 
bell, consistent with the OTP model being less strongly 
correlating. Moreover, it was found that these systems to 
a good approximation have master isomorphs, i.e., that 
all isomorphs have the same shape in the virial/potential 
energy phase diagram. This applies for the total virial, 
but also after subtracting the constraint contribution. A 
general master isomorph was identified between the in- 
vestigated model systems after this subtraction. We do 
not at the present have an explanation for this observa- 
tion. 

A full theoretical understanding of the implications of 
rigid bonds remains to be arrived at. For instance, the 
shape of molecular isomorphs is different from the shape 
of Ref. [l6|for atomic isomorphs. The rigid bonds seem in 
general to increase 7 and decrease the correlation coeffi- 
cient R with respect to the unconstrained system. More 
specific; R decreases significantly with increasing asym- 
metric dumbbell bond length {R w 0.65 around unity 
bond length, see Sec. IIV Al) . This is consistent with the 
results of Chopra et ali^,, who note a worse scaling of the 
reduced relaxation time and diffusion constant with (ex- 
cess) entropy based quantities when increasing the bond 
length of a rigid symmetric LJ dumbbell model. On the 
other hand, it is noteworthy that strong correlation is 
observed for the OTP model even though it has unity 
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bond lengths. The molecular center-of-mass structure in 
reduced units is predicted to be invariant along an iso- 
morph, however, for the OTP model the reduced particle 
structure seems more invariant along an isomorph than 
the reduced molecular center-of-mass structure. The for- 
mer is not predicted to be invariant along an isomorph, 
and thus the observed difference should be investigated 
in more detail to clarify this issue. 
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Appendix A: Constrained dynamics and the virial 



algorithm. The reader is referred to Refs. 113 and [3^ for 
details concerning this aspect. 

The virial W is defined hy W ^ l/3E,=i In 
an atomic system with LJ pair potential interactions the 
virial is given hy W = Wlj = -'^/'^Y^^jriju' {rij). If 
the system has bond constraints -i/)" — {va^i — Va^jY /2 
= ^^/2 = 4 ^/2 it follows from Eq. that the 

constraint force contributes to the virial as Wcon — 
l/3Etir.-G, = l/3EtiA"<,,. 



Appendix B: Constrained NVE and Nose-Hoover 
NVT dynamics in reduced units along an isomorph 

We start our considerations of this section from the 
constrained equations of motion derived in Appendix VK[ 
Eq. (IMl) : 



Constrained dynamics is discussed in many different 
places, see for instance Refs. HI, Hi], and 113 We give 
here a brief introduction to constrained dynamics and the 
connection to the virial expression used in this article. 

Gauss' principle of least constraint^'' states that a clas- 
sical mechanical system of N particles with constraints 
deviates instantaneously in a least possible sense from 
Newton's 2nd law, i.e., that 



mA'Ti - 

i—l 



(Al) 



is a minimum. Here and F, are the position and 
interaction force of particle i. In the case of no con- 
straints, setting the partial derivative d/dvi to zero im- 
plies Vi — Fi/nii — 0, i.e., Newtons's 2nd law. 

In the case of holonomic constraints ^/;"(r^) — where 
a — 1, ...,G, the variation can be carried out by intro- 
ducing Lagrangian multipliers, i.e.. 



N 
i=l 



r F n 2 

L TO,- J 



(A2) 



should be a minimum. Setting the partial derivative 
d/dvi to zero implies (where the factor of one half has 
been absorbed in the Lagrangian multiplier) 



(A3) 



a=l 



Newton's 2nd law thus remains valid if an additional force 
is added (called the constraint force Gi). At this point 
A" is undetermined; however, an explicit expression^ for 
A" can be determined from the condition Vij • icij + r^^ = 
0. In molecular dynamics simulations it is imperative 
to calculate A° correctly to achieve a stable numerical 



mi-ri= Fi 



A^Vr.?/'" =F, + G, 



(Bl) 



Here and F.^ are, respectively, the position and interac- 
tion force of particle i, and A" a Lagrangian multiplier for 
the a'te constraint For simulating rigid molecules^ 
the constraints would in general be a combination of con- 
strained bond lengths = (r^^i — rQj)^/2 = ~ 

c^^ij/^ and linear constraints tp^ = J27=i ^i3i^i — r^j = 0, 
where Cpi is a factor that depends on the geometry of the 
molecule (see Ref. for more details). For simplicity 
we consider only bond constraints in the following. 

A general expression for the Lagrange multiplier A" 
can be derived and is given by^i^ 



N 



N 



/3 = 1 i=l 



1=1 



J ' 

(B2) 
(B3) 



Defining reduced units for length, energy, and mass as 
follows 



U = U/ksT, 
reduced units for time and force follow as 



(B4) 
(B5) 
(B6) 



(B7) 
(B8) 
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Inserting the above definitions in Eqs. (jBll) - (IB3|) and 
using that Vr-V'" = ^a,ij: we arrive at the constrained 
NVE equations of motion in reduced units 

G 

m^-h = F^ + = F, + G„ (B9) 

where 



G N N ^ ~ 

13=1 i.j=l i=l ' 

(BIO) 

N ^ 

Z y^l^^,^2lm_ (Bll) 

Since in general = P^^^c^^ij the reduced constrained 
equations of motion are not invariant along an isomorph. 

Considering instead the molecular center-of-mass mo- 
tion in reduced units 



M, • 'icM.^ = FcM.^, (B12) 



where FcM,i and Mi are, respectively, the reduced force 
and mass of molecule i. Since the reduced force FcM,i is 
invariant along an isomorph, it follows that the molecu- 
lar NVE equations of motion are invariant along an iso- 
morph. The invariance of FcM,i can be seen as follows. 
The isomorph definition Eq. ([5|) implies for a fixed state 
point (1) and arbitrary state point (x), both along the 
same isomorph [where R = {p~^^^rcM.i, Oi, i^i, V'lj 
p^^/^icM,N, On, (I>n, iPn)] 

- U{B}"^)/kBT., = -{/(R^'V^B^i - In Ci,. (B13) 
Taking the gradient V^cJ^^ ^ it follows that 

FS./, = rS/,. (B14) 

This concludes the proof of the isomorph invariance of 
the reduced molecular center-of-mass equations of mo- 
tion. A similar situation is given for the constrained NVT 
equations of motion, however; considering the molecular 
motion, the constraint force disappears, and the proof 
in analogous to the above and shown for atomic sys- 
tems in Ref. [l^ In this case the time-constant of the 
Nose-Hoover algorithm needs to be adjusted along the 
isomorph, otherwise the dynamics is not invariant. 
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